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I. INTRODUCTION 



The development of any new aerodynamic configuration relies heavily on 
Computational Fluid Dynamics (CFD) to enable an economical study of parameter 
sensitivities and configuration integration. The Similarity rule, a more traditional 
technique, was commonly used, prior to CFD, in conjunction with customary methods, 
such as wind tunnel testing and full scale tests. Similarity rule can be applied to cases 
under investigation together with other known information about the test object. 

As part of a design team that was established at the Naval Postgraduate School 
(NPS), the Joint Standoff Weapon (JSOW) missile that is under development is being 
analyzed and geometrically modified to optimize its weight and drag characteristics. The 
aerodynamic results can then be used as background for studying and verifying the Area- 
Rule defined in the 1950’s by Whitcomb [Ref. 1], 

The CFD analysis may be applied using both single-grid and multi-grid techniques. 
State-of-the-art codes (such as GRID GEN Ver. 9.6), can be used to create grids, merge 
them together for physical computational integration and create interpolation regimes, 
while solving flow around the given configuration both in inviscid (Euler) and viscous 
(Navier-Stokes) flows by using CFD codes such as NASA’s OVERFLOW flow solver. 
The multi-grid method allows complicated geometrical configurations to be analyzed; 
however, it is well known that generating grids around configurations that appear simple 
initially, may turn out to be a complicated process. 

The CFD process, apart from involving many primary and secondary codes, is 
known to be very sensitive to the grid generation process and architecture. OVERFLOW, 
for example, is sensitive to the distribution of grid lines. Other sensitivities are found in 
merging codes such as PEGSUS if the interpolation regime boundaries are set improperly, 
resulting in non-overlapped regimes or “orphan” grid points. 
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Only after the CFD process is successfully implemented may the Area-Rule study 
be investigated. This rule, in essence, states that the drag of a body depends (above all) on 
its cross-area distribution along its axis. Pressure related drag decreases as the distribution 
of the cross-section area becomes smoother. 

A combination of Captive Air Training Missile (CATM) and pylon can be used to 
investigate the Area-Rule by running different cases for different CATM-pylon relative 
locations. Methods for achieving area distribution change must be consistent with their 
ease of implementation with the CFD tools, e.g., modeling the area distribution changes by 
grid modification is difficult to implement with CFD; hence, using two grids in different 
axial relative location is much easier. 
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II. BACKGROUND 



A. JSOW WEAPON SYSTEM 

The next generation of air-to-ground (A/G) weapon systems, which are under 
development throughout the world, is a stand-off type of weapon that enables the 
launcher crew to be out of the threatened area, preferably out of ground-to-air missile 
range, while still being able to launch and control the weapon. 

Several such systems are in various stages of development such as SLAM (Stand- 
off Land Attack Missile) developed by McDonnell Douglas (MD) for the U.S. Navy. The 
modified and improved version of the SLAM (SLAM-ER). CASOM (Conventional 
Armed Standoff Missile) developed by the British. JDAM and the AGM 154 (Air-to- 
Ground Missile)-JSOW developed under a joint U.S. Navy/ Air Force contract by Texas 
Instruments (TI) as the prime contractor. 

Three versions of the operational AGM 154 weapon system are being developed 
[Ref. 2]. 



• The basic configuration (AGM- 154 A) cluster weapon contains 145 
Aerojet/Olin BLU-97 submunition (combined effects of shaped charge, 
blast/fragmentation and incendiary) which is designed to replace the anti-soft 
armor Rockeye A/G weapon. 

• The advanced configuration (AGM-154B) cluster weapon contains 6 Textron 
Systems BLU-108. Each contains 4 “Skeets” self-homing, anti-armor 
submunitions. 
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• The unitary configuration (AGM-154C) blast/fragmentation warhead with a 
terminal guidance system. 

An advanced version will include jet engines to increase the range of the weapon 
up to 100 nautical miles. 

Fulghum [Ref. 3] described the weapon as a “truck” that will be able to deliver its 
payload to a nominal range of 40 miles and will be able to be carried by most NATO 
and US aircraft. The future configuration will include smart seekers for the end game. 
The main candidates for the seekers are SAR, Imaging Infrared, Laser Radar, and 
millimeter wave radar [Ref. 3]. The weapon’s operational principal includes the release 
from the airframe platform, gliding to its target using a combination of Global 
Positioning System (GPS)/inertial navigation systems, providing a navigation accuracy of 
up to 10 m, and dispersing the submunitions or hitting the target using the “man in the 
loop” preprogrammed schedule. The development of the weapon system itself is, in 
general, widely embedded with modeling and simulations in order to decrease the 
evaluation tests [Ref. 4], 

B. JSOW TRAINING POD - CATM AND NAVAL POSTGRADUATE 
SCHOOL TEAM INVOLVEMENT 

While the JSOW weapon system is being developed by TI, the Naval Postgraduate 
School (NPS) as a Navy organization has been cooperating with the contractor to develop 
the training configuration of the system. 

The JSOW training program, in addition to having a sophisticated and expensive 
systems such as the SLAM and GBU-15, has degenerated configurations that are being 
used for training purposes only. The CATM configuration is a pod that hangs on the 
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regular aircraft pylon, simulating the flight condition of the operational version, and has 
characteristics similar to the operational vehicles. 

A working group consisting of professors and students from the Aeronautical and 
Astronautical Department at NPS was established in 1995 to develop a training pod 
including aspects such as: 

• Software simulation of the JSOW training mission. 

• Static and dynamic load analysis. 

• Weight optimization (to achieve the 300 lb limit that obviates the need for 
emergency separation). 

• Quantity in fleet optimization for the training pod. 

• Aerodynamic analysis for the pod. 

• Drag optimization. 

The group started a preliminary development of a training pod that would achieve 
the 300 lb limit and have a relative low drag, thus, enabling longer training missions. 

This thesis concentrates on the aerodynamic analysis aspect of the training pod, 
particularly on the drag optimization of the CATM. 

C. CFD - COMPUTATIONAL FLUID DYNAMICS AS A DEVELOPMENT 
TOOL USED IN THE PROCESS 

1. Overview 

CFD is a powerful tool for aerodynamic analysis. Its use has accelerated 
significantly during the last few years as a result of the fast development of 
supercomputers, which produce solutions within a reasonable time and effort. This by the 
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solution to an enormous number of simultaneous differential equations which are needed 
to analyze flow around complicated vehicles. 

2. Competitive Methods 

Competitive methods for aerodynamic analysis such as wind tunnels, pure 
analytical and parametric applications, and flight tests suffer from high cost limitations, 
inaccuracies, and availability difficulties. Ref. 5 quotes the following: 

The Navy requires a predicted method capable of predicting the 
store’s pitching and yawing moments to within 10% of flight data in the 
transonic flight regime. 

Assuming CFD is accurate in store’s pitching and yawing moments calculations 
to within 10% of flight data in the transonic flight regime, it should satisfy the 
requirement above. 

D. THEORETICAL BACKGROUND 
1. General 

OVERFLOW, the flow solver, developed by NASA Ames was used in this 
research, it is capable of solving the Euler as well as with Navier-Stokes (NS) equations. 
These are the inviscid/viscous representations of flow characteristics for any 
configuration defined by initial and boundary conditions. The detailed characteristics of 
the code will be discussed later. There are no closed solutions for Euler/Navier-Stokes 
equations. The analytical closed solutions are restricted to the potential flow theory, e.g., 
the solution of a flow past a sphere or a flow past a cylinder. 

In essence, OVERFLOW is oriented to solve fully viscous flows. The Euler 
solution is a degenerated form of the full solution achieved by turning off some options 



6 



from the full NS solution. The viscous terms, partially defined by the shear expressions 
are set to zero in an Euler case. 



The code is an implicit, conservative solver that can use the Chimera overlapped 
grid scheme in conjunction with PEGSUS. 

2. Navier-Stokes (NS) Equations And Solution Techniques 

The NS equations include the viscous effects concentrated mainly in the boundary 
layer [Ref. 6]. Thus we can introduce the NS equations in their Cartesian form: 

Introducing the following variables: 

• u, v, w - the flow velocities in the x, y, z directions respectively. 

• p - the flow density. 

• p - the pressure. 

• T - the shear stress tensor. 

• e - the internal energy. 

• E t - the total energy. 

• q - the heat flux vector. 

• k - thermal diffusivity of the fluid. 

• p - viscosity coefficient. 

• T - temperature. 

• R - universal gas constant. 

• C v - specific heat at constant volume. 



BQ BF BG BH 
— + — + — + — 
Bt Bx By Bz 



0 ) 
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Where Q, F, G, H are defined as follows: 



Q = 



P A 
pu 

pv 
p w 



( 2 ) 



pu 

pu 2 +p-x„ 

F= puv-t xy 

P UW-T „ 

(E, + p)u-ux xx - VX xy - WT„ + q x J 



(3) 



G = 



pv 

PUV-^y 

pv'+p-Tyy 

p vw-x ^ 

(£,+p)V- UT.y-VTyy-WT^+qy^ 



(4) 



H = 



p w 

puw-x M 
pvw-x^ 
p w 2 + p-x zz 

(E, + p)w-ux xz -vx yz -wx zz +q z , 



(5) 



Where the shear and heat flux components are defined as: 



x 
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2 ( d u d v 

3 ( 3x d y 
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(7) 

( 8 ) 

(9) 

( 10 ) 
( 11 ) 
( 12 ) 

(13) 

(14) 



And the equivalent equations of state, internal energy and total energy are: 



p = p RT 



e = C v T 



(15) 

(16) 



E t =P 



e+ — (u 2 + v 2 + w 2 ) 
2 



(17) 
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m. COMPUTATIONAL TOOLS AND TECHNIQUES 



A. ARCHITECTURE 

One of the main characteristics of CFD work is the complexity of various codes 
involved in the analysis. It was discovered that certain aspects of CFD, especially grid 
generation, tend to require more “art” than science to achieve reasonable results. 

In this section we will describe the codes that were used in the process and the 
relationship between them. The techniques for using them is described in detail in 
Appendix C and Appendix J. 

List of codes (not necessarily in the order of use): 

1. GREDGEN - Grid generation program 

2. PEGSUS - Determine the interpolation stencil (IB LANK file) 
between multiple overlapping grids. 

3. OVERFLOW - Row solver program. 

4. FAST - Visualization program of flow solution over the 
investigated blocks. 

5. GRIDED - Editing grids program. 

6. RMG2PEG - Merges the separated grids into one input file for 
PEGSUS. 
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7. MERGE41 - Merges the output files of PEGSUS and IBLANK file 
into the “grid.in” input file for OVERFLOW. 

8. xa2b.f - ASCII to Binary and “xb2a.f \ Binary to ASCII transform 
files. 

9. readfomo.m and readresid.m - Matlab files (written by Bret S. 
Barton [Ref. 7] and modified in this research) makes it possible to 
present graphically the OVERFLOW output data given in 
“fomo.out” and “resid.out”, respectively. 

10. wrq - a conversion file that converts the solution file “q.save” from 
its Cray format to a SGI format. 

11. wrgr - a conversion file that converts the grid file - “grid.in” in its 
Cray format to a SGI format. 

B. GRIDGEN 

1. General 

The basic element of the CFD analysis is modeling the body under investigation in 
its environment. Currently, NPS uses the following programs to create meshes: 

a. GRIDGEN 

b. HYPGEN 

c. GRAPE 
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The code that was used in this research was GRIDGEN, developed by John P. 
Steinbrenner and John R. Chawner, MDA Engineering, Inc. [Ref. 8]. Version 9.6 of the 
code was used throughout this research. Currently, version No. 1 1 is under development. 
The original version of the code was developed by the U.S. Air Force. Since then, the 
code became commercial, and attempts by NPS to update that version (instead of the 
current 9.6 version, which contains some bugs) have not been successful. The code was 
developed using the Silicon Graphics IrisGL graphics library [Ref. 8], and hence can be 
used by the computers in the computer laboratory of the NPS Aeronautics and 
Astronautics Department. The code is menu driven, very user friendly and offers a very 
convenient way to create an optimal grid around complicated surfaces. Difficulties that 
were found in the code will be elaborated further on. 

2. Theory 

In essence, the purpose of the code is to create a block of a large amount of 
volumetric cells that fill the investigated domain. 

Two domains considered in the grid generation and flow solver processes: 

• Physical domain, in the grid generation program. 

• Computational domain, in the flow solver program. 

The physical domain (assigned by x, y, z coordinates) consists of the actual 
surfaces that can be modeled geometrically. The computational domain (^,r|,Q is a 
virtual domain of simple cell boxes that can be used easily by the flow solver [Fig. 1]. 
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Figure 1. Physical and computational domains used by GRIDGEN [Ref. 6]. 

The output of the grid generation process is a file that contains all the coordinates 
of the intersection points of the grid lines. 

A grid (or mesh, as it is called in other sources) is a block made of six domains (as 
a box). The domains are made of four connectors. The basic element is a connector, that 
is, a line that can form any desired shape and is divided into a certain number of grid 
points. A connector [Fig. 1] can be made of single or multiple segments. This architecture 
provides a convenient way to model complicated shapes by dividing the shape contour 
into many short, accurate segments. The connector is the subelement of the next grid 
element, the domain [Fig. 1]. The domain is a surface made of 4 connectors: 2 pairs that 
create a 3D mesh of the surface. In order to achieve a reasonable domain, each two 
opposite connectors must be the same dimensions (but not necessarily the same grid 
distribution). In the domain menu options [Fig. 2], one can identify the power of 
GRIDGEN to create complicated, physical domains, e.g., by default, GRIDGEN provides 
algebric solver to smooth the grids; however, in complicated configurations, the feature of 
the Elliptical Solver makes it possible to curve grid lines in the appropriate direction and 
create a uniform mesh, thus preventing negative cells in OVERFLOW (which will be 
elaborated in the following sections). 
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Figure 2. Domain main menu in GRIDGEN. 



The domain is the subelement of the upper level element of GRIDGEN, the block. 
The block is in essence a box made of six facets: the six domains. To achieve a valid 
block, the dimension of each two opposite domains has to be identical [Fig. 1]. GRIDGEN 
makes it possible to create a block with more than six domains so that the final 
configuration after summation (joining) of the appropriate domains will be valid, i.e., the 
grid surface dimensions of each two opposite domains will be the same. In all the blocks 
that were created in this research, the configuration included one domain that was the 
investigated domain (i.e., the surface of the CATM or the surface of the pylon, or the 
surface of the wing became one domain in each of their respective blocks). The other five 
domains were built around them to create a block. 
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3. 



Difficulties 



a. General 

One of the most problematic phases of this work was the task of grid 
generation. As was mentioned above, some aspects of CFD can be identified as an “art”, 
and grid generation is at the top of this group. 

The difficulties can be divided into three main categories: 

• Internal inherent bugs, 

• Difficulties caused by complicated geometry, 

• System problems. 

Concerning the first category, we have been in communication with the 
GRJDGEN developer, and have described the problems and irregularities that were 
identified during this research. Since the current version under development is No. 1 1 and 
no technical support was provided for the NPS version ( No. 9.6), we had to reconcile the 
two versions and find ways to overcome the obstacles, as will be described below. 

b. Internal Inherent Bugs 

(1) Grid collapsing. Most of the grids used in this research were of 
the “O” type grid, which is usually used for symmetric objects. All “O” grid connectors 
are located on one surface. In a few circumstances, the grid/block collapsed to a zero 
volume configuration and the work had to be redone. 

This phenomenon occurred in several cases: the first case 

happened while attempting to save a block after it was created by adding six domains 
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together. The “save” command caused the grid to collapse. The solution to this problem is 
to create the six domains first, save them as a GRIDGEN/ ASCII file, reload this file, and 
create the block. No logical reason was found for either the problem or the solution. 

The second case happened while attempting to modify the grid 
(e.g., to redimension or redistribute a connector). This caused the grid to collapse as well. 
No solution was found in this case, except to create the block again after modifying the 
connectors/domains as connectors/domains files (for that reason a copy of the connector 
and domain files must always be saved). 

(2) Joining domains. The basic requirement for joining two grids in 
general, and particularly two domains together, is to maintain the grid dimension 
compatibility of the domain connectors being joined. However, in a few cases it was 
found that, although the requirement was fulfilled, the code refused to join or agreed to 
join the unwanted connectors (e.g., in case of multioptions of joining configurations). The 
solution for this is to create one of the two domains again. This time, the connectors 
should be added in the direction opposite to the original order, i.e., if the original domain 
was created by adding the connectors in clockwise direction, it will be created by adding 
the connectors in the counter-clockwise direction. Otherwise create one of the domains 
again, which, this time, should be started with another connector than was originally 
started. 



It must be mentioned that it is written specifically in the manual, 
and the developer insists that the domain generation is insensitive to the order of adding 
connectors, but, as detailed above, the results of this research were otherwise. 

(3) Creating blocks. There are cases where an attempt to create a 
block fails and a “wrong topology” error message is given, even though the block is made 
of perfectly compatible domains in dimension aspect. In this case, it is recommended to 
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check the direction of £,r|,£ in each opposite domain, even though, the domain 
characteristics are not supposed to be sensitive to the direction of creation. If the direction 
of one of is opposite to the direction of the appropriate in the opposite domain, the 

code might refuse to create the block due to this reason and one of the domains in the 
couple has to be rebuilt in the opposite direction. 

c. Difficulties Caused By The Complicated Geometry 

The fewer the grid blocks involved in the analysis, the less complicated the 
analysis, but, on the other hand, the less accurate. 

Since the major configuration investigated was a combination of a missile 
and a pylon, an attempt was made to create one block that would model the missile and 
the pylon as well. The main reason was to prevent, at least in the beginning of the 
research, the need for using PEGSUS which contributes a significant amount of 
complexity, as will be described later. 

A single block grid could not be successfully generated because of the 
complicated configuration, mainly because of the difference in the length of the missile 
compared to the pylon, and to the curvature contour of the pylon, which caused a very 
dense grid in a certain region of the pylon [Fig. 3]. 
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Figure 3. A dense grid structure of a joined missile/pylon block. 
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The phenomenon that occurred due to the very dense grid regime is the 
“Negative Volume” cells computed by OVERFLOW. When the flow solver identified 
negative cell volume, once it had read in the grid, the flow solution was not advanced and 
the code stopped with an error message. Using the smoothing tool of the elliptical solver 
and a very coarse grid finally produced a grid that was acceptable by OVERFLOW, but 
was physically worthless [Fig. 4], 

Fig. 3 presents another attempt to get rid of the negative cells. In this case, 
the investigated domain was copied and the copy was placed next to the investigated 
domain in a way to assure that all the grid lines are parallel and, hence, no negative cells 
can be created. Although, theoretically, negative cells should not be in this configuration, 
they appeared. 



d. System Problems and Difficulties 

(1) Segmentation fault. The main system limitation of GRDDGEN 
processing is “Segmentation Fault” errors which are caused by the system and probably 
have no direct connection to GRIDGEN. They are caused by the network’s inability to 
handle heavy files. These phenomena occur when trying to load a certain file (grid), or 
even when doing a simple action such as domain creation. The code is killed and a 
“segmentation error” is given. Usually it takes a while to get rid of these segmentation 
faults, since they are a problem with the network rather than the code. 

(2) Problems caused due to file storage in the Cray. Since most of 
the research was done using the Cray, and due to the large grid and solution files (order of 
20M for a file), the files were stored on the Cray disks and tapes. To efficient the capacity 
management of the files on the Cray, the files are migrated after a certain time to storage 
tapes (usually after 48-72 hours if the files were not used). It was found that in some 
cases, due to network difficulties, files that were stored on disks, could not be read 
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properly from the Cray by simply changing directory to the Cray from the SGI machines. 
Those files were loaded as corrupted files, embedded with symbols not written in the files 
originally. In those cases, the only solution is to transfer the files to the local machine 
(SGI) by file transfer procedure and then use them. Another difficulty caused by this 
“migration” method is that the time to load a migrated grid or solution file is very long, 
and can take 10-20 minutes. This phenomenon causes delays in the analysis process. 
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Figure 4. A coarse grid structure of a joined missile/pylon block. 
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c. 



PEGSUS 



1. General 

The PEGSUS code enables the investigation of complicated configurations 
containing as many blocks as needed, hence preventing the need to create complicated 
structures that can be divided into simple ones. 

The current version being used at the Aeronautics computer lab is 4.0 [Ref. 9], In 
this research, the PEGSUS code was used in conjunction with the flow solver, 
OVERFLOW, through intermediate code, MERGE41, and with the Chimera method, to 
interpolate between each two adjacent grids. 

2. Chimera Method 

The concept behind the Chimera method is given in Ref. 9 and described in Figs 5 

and 6. 




pylon mesh 



Figure 5. PEGSUS / Chimera concept - mesh to mesh communication [Ref. 9]. 
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Interpolation boundary 
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Interpolation boundary 
point in pylon mesh 




Hole creation 
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^-Hole boundary 
n CATM mesh 



Hole in CATM 
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Figure 6. PEGSUS / Chimera concept - overlap region between meshes [Ref. 9]. 

Here are the main characteristics of the Chimera technique: 

Each mesh of the two mesh couple receives and donates flow field information to 
and from the other. The interaction between the two meshes takes place in an 
“Interpolation Regime”, which is embedded in both meshes. 

To create this interpolation regime, two boundaries must be defined, which creates 
three different regimes: 

• Hole point regimes where the donated grid points have to be computed, but not 
the accepted grid points. 

• Interpolation regime where the information has been transferred from one grid 
to another. 

• Field regime where the accepted grid points have been calculated with no 
connection to the donated information. 
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3. Technique 



Once the individual grids are available from GRIDGEN, the merging of the grids 
can be accomplished. As a preprocedure to PEGSUS, one has to run the RMG2PEG code 
(see details in Section G ) to create an input file to PEGSUS: “INGRID”. 

An input file is created [Appendix B], which in principle, contains descriptions of 
the grids, including which grids are linked together, which range of grid points has to be 
included in the updated regime ( J, K, L INCLUDE variable), and a detailed definition of 
the HOLE boundary. What grid the HOLE is part of (ISPARTOF) and in what grid the 
HOLE is made in (MHOLEIN variable), as well as a precise definition of the surfaces. 

4. Orphan Points 

While the main obstacle in OVERFLOW/GRIDGEN is negative volumes, the 
major problem found in PEGSUS is Orphan Points (OPs). 

The simple definition of an OPs as given in the user’s manual [Ref. 9] is, “Any 
boundary points that remain after the mesh priority list is exhausted are termed 
‘orphans’”. The manual also mentions that, in the flow calculation point of view, an OPs 
is treated as a hole point, i.e., it does not get updated by the flow solver. 

If the number of OPs is limited, one can ignore them and “pay the price” of a few 
unupdated grid points; however, the number of OPs is usually more significant and they 
must be eliminated or significantly reduced (see sensitivity study in the following 
sections). 
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Figure 7. Orphan Points definition [Ref. 9]. 

To eliminate the OPs, few methods are recommended. 

• Grid shifting. 

As described by Barton in Ref. 7, a slight shift between the grids can significantly 
eliminate the number of OPs; however, this procedure might be very tedious. Since it is a 
trial-and-error procedure in each iteration, one of the grids has to be modified, as 
described above. Modification of a grid (especially using GRJDGEN) can be very time 
consuming. 
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Enlarging the interpolation regime 



• 

In the PEGSUS input file, the inner and outer boundaries of the interpolation 
regime are specifically defined. In general, as the boundary increases in size, the 
probability of getting OPs is lower since the probability that a point won’t have any 
adjacent point to get information from is low. 

• Decreasing the range of included calculated grid points 

Quite often, the OPs appear in the outer layers of the block. Since the outer layers 
are in the far field, one can limit the calculations to the whole grid, except for a few 
(approximately 2 to 3) grids on the outer boundary. It should be mentioned that those grid 
points are not updated by the flow solver and one has to check carefully if it does not 
effect the whole solution. 

A sensitivity study has been done to investigate how the flow pattern can be 
influenced by the existence of OPs, the results of which are detailed in the multi-grid 
section ahead. 



5. PEGSUS Output 

PEGSUS creates several interpolation files, for example, “fort.2”, as well as data 
files such as the OPs list and the EBLANK array (which contains “0” for hole points and 
“1” for field points). 

The main file is the “grid.in” that, with the interpolation data are the input files to 
OVERFLOW the “grid.in” file is created actually after MERGE41 is applied on the 
PEGSUS output files. 
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D. OVERFLOW 



1. General 



The FORTRAN code, OVERFLOW, was developed as a flow solver code capable 
of using either the full Navier-Stokes or only the Euler description. Using the Chimera 
method in cooperation with PEGSUS and GRIDGEN, the code makes it possible to deal 
with complicated configurations made of several separated blocks. 

At NPS, Version 1.6aw is installed on the Cray as well as on the SGI machines. 
The code was run in both configurations, as shown in the following sections. 



2. Technique 

The input file [Appendix B] includes the variables needed for the code. The 
detailed characteristics of the code are given in the OVERFLOW user’s manual [Ref. 
10]. Here are some highlights: 



The principal flow data (e.g., Mach number, AOA, Reynolds number, temperature 
at infinity) are given in the “FLOINP” card. The next cards include information 
concerning the smoothing, techniques of integration, time steps, and so on. 

Two important values are given in the “TIMACU” card: DT, the time step of the 
integration, and the CFL (Courant-Friedrichs-Lewy number). 



The CFL is given by the following expression [Ref. 11]: 

ax At 



N = 



Ax 



(18) 



where “a ” is the wave speed of the following first order equation: 
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(19) 



which is the degenerate form of the Euler basic equation without the source terms. 

The convergence and stability of the solution depends, to a large extent, on the 
combination of these two values, DT and CFL. One of the major stages in the analysis is 
to fmd the optimized combination, which may be different from one case to another and 
depends on the Mach number as well. 

The final section of the input file deals with the boundary conditions (BCs) to be 
applied. There are different BCs for viscous and inviscid solutions. The BCs can be 
defined specifically on certain surfaces using the initial and ending grid indices of the 
desired regime or on the whole grid using [-1,1] as min / max. definitions. 

After the input file and the appropriate “grid.in” file (including interpolation files 
in the case of multi-grid) have been created, the code is run and generates several output 
files: “fomo.out”, a history file of all the aerodynamic coefficients calculated in each 
iteration, “resid.out”, residuals calculations (for each iteration) and solution file, “q.save” 
(contains information of the “Q” vector defined in equation 2 for each grid point). 

In the case of a non-convergent run or a numerical problem that causes a negative 
calculation of the pressure or the density, a “q.bomb” file will be generated. This file 
contains the last valid information of the “Q” vector and can be analyzed as “q.save” 
being. Finally, more information concerning the values of the pressure and density is 
given in the “rpmin.ouf ’ file. 

The “q.save” file can be visualized using FAST and the tabulated information in 
“fomo.out”, “resid.out” files can be plotted using Matlab codes READFOMO.M and 
READRESID.M written by Barton [Ref. 7], 
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3. 



Useful Suggestions 



a. Unexplainable Negative Cells 

Cases where a negative cell error was given occurred quite often even in 
“restart” runs with the same input file and grid that was run the first time when the error 
was not given. The best solution is to copy the “grid.in” file again to the directory and run 
the code again. 

b. Viscous Results While Setting The Viscous Terms To “F” 

One must keep a consistent connection between the BCs codes and the 
flow definition, i.e., for example, if one sets the viscous terms to “F’ ( in the VISCJ, K, L 
parameter), and maintains the boundary condition code as a viscous one (“5”, the viscous 
adiabatic wall), the friction data presented in the “fomo.out” file will not be zero as 
expected in the inviscid flow, but it will contain residual small values (2 to 3 orders of 
magnitude less than the pressure drag). 

E. FAST 

The visualization program for OVERFLOW solutions given in the “q.save” file is 
FAST (Flow Analysis Software Toolkit). This code, when used with powerful SGI 
machines, gives an impressive 3D way of visualizing flow characteristics. The code 
receives the grid files and the solution file as an input, and calculates the value of the 
desired parameter in each cell. 

In case of a multi-grid run, the “grid.in” file created by the MERGE41 code must 
be used as the grid file (in case of Cray run, after converting the “grid.in” to a file 
readable by the SGI machines by using “wrgr” code) and not the original separated grid 
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files because, in the latter, the blank grids won’t be seen and a non-physical solution will 
be presented. It should be stressed that FAST will present the solution file “q.save” on a 
grid file whether it contains blank points or not. Hence one must specify to read the blank 
points information (in the bottom of the grid file). 

A few useful suggestions should be mentioned concerning the code: First, since 
the code runs on the SGI machines, it is very comfortable to run the initial, non-CPU 
consuming cases on the SGI machines and check them directly by FAST, i.e., “q.save” 
generated by OVERFLOW, which runs on the SGI can be read directly by FAST (as a 
“solution”, “unformatted” file). On the other hand, in the Cray case, the “q.save” file has 
to be converted to a SGI readable format using a “wrq” conversion file (see next 
sections). Running “wrq” can take a while. Similarly, the “grid. in” generated by the SGI 
“xa2b” file can be read directly without any conversion, again unlike the Cray “grid.in” 
file that has to be converted by “wrgr” code. 

Another important point is that FAST does not care about the grid distribution, 
i.e., one solution file can be imposed on two grids that are identical in the grid dimension 
but completely different in the grid distribution. This might cause confusion, unless 
careful follow up is being done on the files. 

F. GRIDED 

A very useful tool for editing grid files is the GRIDED program. This menu- 
oriented file enables one to change the orientation of the grid direction, flip between 
indices and more. This is the code used mainly to edit new blocks to a uniform 
configuration (in the I, J, K directions and orientation point of view) of new blocks. It’s 
importance lies in the fact that, once the orientations and dimensions are the same, the 
same input files can be used for OVERFLOW. 
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It should be mentioned that this version of the code is probably not the latest one. 
In the “test” directory of OVERFLOW where test cases are given, few of the grids need to 
be “expanded” using GRIDED, for example, the axisymmetric test case; however, the 
expansion must use option ”8” in GRIDED which is absent in the NPS version, and hence 
not all the test cases could be run. 

G. RMG2PEG 

This intermediate file creates the merged grid file for PEGSUS, it is a menu 
formatted file and uses the unformatted grid files [Appendix B]. 

The file simply adds together the separated grid files one on top of the other and 
creates one unformatted grid file called “INGRID”. 

H. MERGE41 

MERGE41 is second intermediate file that merges the output files of PEGSUS 
into a PLOT3D format to be read by OVERFLOW. 

The file uses the blanking scheme produced by PEGSUS and gives a “grid.in” file, 
which contains all the grid to be run in OVERFLOW including the blank grid points. 

The hole can be seen well in Fig. 8, where the pylon makes a hole in the CATM 
grid, which can be identified by the hole made in the CATM grid through which the pylon 
can be seen. 
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Figure 8. Hole creation in the CATM by the Pylon. 



I. 



XB2A.F, XA2B.F ■ BINARY TO ASCII AND ASCII TO BINARY 
CONVERSION FILES 



Since the flow analysis process uses formatted and unformatted files in different 
applications, these two conversion files change the format of the desired grid. First one 
has to rename the grid being investigated to “grid.for” (when running “xa2b”) and get the 
unformatted file - “grid.in”, notice that if “grid.in” file is already exist in the directory, the 
code will fail - make sure to erase the previous “grid.in” file. Similarly, in the “xb2a” file, 
one has to rename the unformatted file to “grid.in”, and retrieve the “grid.fmt” formatted 
file. 



As mentioned above, the code was installed both on the Cray and the SGI 
machines. In order to run the code, an unformatted grid file had to be created. The file 
was generated by the conversion file “xa2b.f” (Appendix B). In general, FORTRAN 
compilation on the Cray will create an executable file, “xa2b”, which will not run on the 
SGI machines and vice versa. Compilation on the SGI will create an executable file for 
the SGI only. 

J. READFOMO.M, READRESID.M 

These two Matlab files were written by Barton for his thesis project [Ref. 7] and 
were modified in this research. The files read “fomo.out” and “resid.out” files which were 
renamed to “fomo_l.out”, “fomo_2.out”, etc. and “resid_l.out”, “resid_2.out” 
respectively in “restart” runs. The original files could read up to four 
“fomo.out”/“resid.out” files. The files were then modified to read up to ten output files, 
and to compare between several cases (e.g., viscous and inviscid runs) 

These files present a very convenient tool for graphing or processing the 
OVERFLOW output files; however, one can use easily GNUPLOT or XMGR, which are 
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also installed in the SGI machines to present the results as well. In this research, results 
were plotted by the Matlab files using the script file “results”, which contained the 
procedure for creating the Cd and residuals as a function of the number of iterations and 
the OVERFLOW input file to have a convenient follow up of the results as a function of 
the input data. 

K. WRQ 

This file, written by Prof. Ismail H. Tuncer from NPS, converts “q.save” from its 
Cray format to a SGI readable file. It is a menu-oriented, robust conversion file that 
enables one to convert a single block as well as multiblock solution files by defining how 
many blocks (grids) are included in the solution file. The file creates another file, “q.cgi”, 
which can be presented by FAST. The code, in its current format must be run on the Cray, 
otherwise, if has to be compiled again. 

L. WRGR 

Similar to the previous conversion file, “wrq”, this file, written by Prof. Ismail H. 
Tuncer from NPS, converts the “grid.in” file in its unformatted form to a SGI readable 
file called “gr.cgi”. This file can be loaded into FAST instead of loading all the sub grids 
separately. This file has it possible to visualize the grid files as they were manipulated by 
MERGE41, and to include the IB LANK points. 
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IV. TEST CASE CONFIGURATION 



A. GENERAL 

In order to study the OVERFLOW features and setup files, few test cases were 
modeled. The models description given in Appendix A, including a description of the 
model and the reasons for developing it. 

B. THE MODELS 

The first test case model (test_4_100.grd) was a symmetric model made of two 
rounded slices [Fig. 9]. Few cases were run using this grid, but the solutions were not 
satisfactory. The assumption was that the sharp comers created by the investigated 
domain with the symmetry upper and lower planes may have caused difficulties with the 
symmetry boundary conditions. The reason for it is that in case of symmetry boundary 
conditions, the flow information on the symmetry plane is being evaluated according to 
the two adjacent planes (one from each side of the symmetry plane), and once the comer 
is too sharp, difficulties in the interpolation between the two adjacent planes are being 
found in the flow solver. For that reason, an axisymmetric test case was modeled [Fig. 
10], and the solutions improved significantly. The grid distribution was carefully done, 
since it was found that using a grid that is too coarse near the surface will cause a 
negative pressure drag coefficient (Cd). Hence it is important to make the grid 
distribution as dense as possible in the boundary layer regime. On the other hand, a too 
dense grid caused negative cells. For example, GRIDGEN/OVERFLOW could not handle 
the grid (bb_dense_4_100.grd) that its critical grid distribution on the “stings” domains 
was as follow - AS initial =0.01, AS final=150 even with 10 iterations of elliptical solver, 
the same grid with AS initial=0.05, AS final=100 and two elliptical solver iterations was 
successful. 
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Figure 9. Test case - symmetric made from two slices. 







Figure 10. Modified test case - axisymmetric configuration. 
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c. 



OVERFLOW 



As mentioned above, the test cases were developed to study the OVERFLOW 
features, and few principals ideas were found that have to be mentioned: 

1. CFL, DT And Convergence 

The values of DT and CFL are critical parameters to the convergence process of 
the code. The sensitivity to the values of the CFL and DT was carefully studied, and the 
OVERFLOW user’s manual [Ref. 10] provides some recommendations concerning the 
initial values for CFL and DT, which depend on the Mach number regime. It was found 
that, as mentioned in the manual, a small value of DT must be used for the first trials, at 
least for the first few tens of iterations. Then, if the solutions seems to converge, the DT 
can be increased gradually. 

2. GRID Topology 

a. Grid Points Distribution 

Besides the importance of a dense grid near the body surface, it was found 
that an insufficiently dense grid can cause such effects as negative Cd, and a lack of 
convergence, besides non-physical solutions. 

b. Symmetric Blocks 

When using an “O” type grid, one must not forget to build the grid in a 
geometry that makes it possible to use symmetry BC; i.e., in the case of symmetric BC, 
there must be a source of information for the surface on which the BC is implied. 
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Therefore, the grid must continue beyond the plane of symmetry, in other words two extra 
planes must be added on the other side of the symmetry plane (see Fig. 11). 

D. CONVERGENCE 

Convergence of the OVERFLOW solution, unfortunately, does not imply a valid 
solution. In many cases, such as the negative Cd cases, the solution converged rapidly, but 
the results were wrong. 

The main parameters that control the convergence are the time step and the CFL. 
Running the first iterations with small values of DT (approximately 0.001 to 0.005) is 
recommended while leaving the CFLMIN as 0.0 (default). The OVERFLOW manual 
recommends [Ref. 10, p. 34] that CFLMIN=5 for subsonic until it becomes low 
supersonic, and that the starting value of DT=1 sec for 50 iterations, which will be 
decreased by a factor of 10 until the solution ceased to “explode”. One must look at the 
behavior of the residuals in order to see if it tends to grow or not. 

As the time step decreases, with no connection to the convergence, the AQ value 
decreasing. (The 8th parameter in “resid.out” represents the residuals, which are different 
in the “Q” vector from one parameter to the next). Hence, one can be confused and get 
“good” values of AQ (approximately 10E-07 to 10E-08) if small values of DT are implied 
whereas the solution really did not yet converge. 
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Upper extra plane (green) 



Lower extra plane (green) 




Upper symmetry plane (red) 



Axi- symmetric surface 



Lower symmetry plane (red) 



Figure 1 1. Two extra planes added to achieve symmetry boundary condition. 
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V. JSOW/CATM BODY CONFIGURATION 



A. GENERAL 

The basic configuration in the analysis is the missile itself, without an additional 
pylon or wing (Fig. 12). The analysis was begun by using the geometry provided by the 
contractor (TI). The geometry was given at first as a set of drawings, since the accuracy of 
the downloading of measurements from the draws is limited. A numerical data set was 
generated and a cross sectional file was sent from TI. The file is in IGES format, and 
attempts to convert this format to a readable format by GRJDDGEN were not successful. 

The grid cross-sections can be loaded as an IGES files but cannot be modified, 
i.e., in order to get a suitable grid one has to create a block around the basic missile 
configuration. This was found to be impossible even by intermediate codes such as 
IDEAS or AUTOCAD. In a parallel research that was done on the structural aspects of 
the CATM, the file could be read into IDEAS without much effort [Ref. 12]. 

An “O” type grid was created around the missile that was modified during this 
research (see Appendix A for details about the various blocks). The first configuration 
was similar to the operational configuration (160 in long) but without its wings and 
control fins. 

After an optimization to reduce weight below 300 lb (Ref. 12) a shorter 
configuration was modeled. In essence, this configuration is the original configuration 
with a portion removed from the constant cross-section portion of the missile [Fig. 13]. 
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Figure 12. Original missile grid. 
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Figure 13. Modified CATM grid after was shortened to reduce drag. 
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B. GRID GENERATION 



1. Methodology 

In general, grid generation has to be made in such a way as to evaluate the flow 
around it, including all associated phenomena. It must also prevent the block boundaries 
from influencing the flow characteristics. All this, given the necessity to “save” as much 
as possible in the computational parts that tend to be high CPU consumers. 

Assuming the flow around the body is symmetric across a central vertical plane 
passing across the body, through most of the research the grid architecture represented 
half the configuration. This assumption is not completely valid, since the actual CATM is 
mounted on an airplane station that has certain flow characteristics that probably can not 
supply a complete symmetric flow around the CATM (due to adjacent stores, aircraft 
structure etc.). This approach is well accepted in CFD and applied to the modeling of a 
full aircraft as well [Ref. 6). 

The first blocks generated were relatively small compared to the dimension of the 
missile, i.e., the diameter of the outer domain was approximately twice the CATM length. 
These configurations caused wave reflections as described further on. Later, the relative 
block size was modified. 

2. Grid Lines Distribution 

In general, it is obvious that the best CFD results are achieved with a smaller 
individual grid cell size to achieve smoothness in the flow calculations; however, since 
the computations use a lot of computer resources, the dense grid orientation is important. 



46 



Dense grids have to be designed in regimes of rapid changes [Ref. 6], e.g., in the 
CATM case, the regime near the body boundary where the boundary layer influence is 
major (Fig. 14). In other cases they must be designed near sharp comers regime (not 
found in CATM case). 
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Figure 14. Grid structure of the CATM - zoom in to the boundary regime. 
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Another situation where a dense grid has to be applied is to create uniformity in 
the block. In order to prevent phenomena such as grid surfaces bending, the grids can be 
redistributed on the main domains. The advantage of this, besides smoothing the grid 
shape, is preventing negative cells in the grid that are caused by extremely bent gridlines. 

Technically, the grid distribution is controlled by the “Redistribution” function in 
the “Connectors” menu. In principle, the code makes it possible to spread the grid points 
along the connector by several schemes, in our case we used the geometrical distribution 
which distributes the grid points by definition of the initial distance (AS) and the final 
distance between them, the gaps lengths are distributed geometrically. It should be 
mentioned that, if the initial definition of AS is “0” - the Jacobean will be zero since the 
volume of the first grids cells is zero. However, GRIDGEN prevents this and sets a 
certain small number (but usually not small enough) to get a valid solution. Hence, one 
must specifically define the initial value of AS and make sure it is not zero. 

C. FLOW ANALYSIS 

1. General 

Getting a solution in a single grid case is relatively simple, but the validity of the 
solution has to be carefully assessed. The CATM model was run in a Mach range of 0.5 - 
1.8. As the Mach number became higher, the solution converged more rapidly due to the 
fact that the influence regime of the flow field became more limited. The Mach angles in 
supersonic flows are small and the interactions with the block boundaries are almost 
invisible (Fig. 15, 16 and 17). In the case of the lower Mach numbers, the convergence is 
slower but can still achieved. 
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In general, the CATM body investigation was the base study for more complicated 
configurations. Therefore, a wide variety of flow conditions were applied such as 
changing the Mach number, angle of attack, and solution type (Euler, vs. NS). 

2. Boundary Conditions 

Overflow offers a variety of BCs for Euler inviscid analysis as well as Navier- 
Stokes solutions. A detailed table containing all the available BCs is given in the 
OVERFLOW user’s manual [Ref. 10], Table 3.3. 

The BCs that were used in this research are as follows (see Fig. 18): 
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Figure 15. CATM flow field, free flight Mach 1.2. 
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Mach, invicicl 



0.534667 



0.754980 



0:975293 



Figure 16. CATM flow field, free flight Mach 0.85, AOA=0 deg, inviscid case 
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0.000000 0.513858 1.027716 



Figure 17. CATM flow field, free flight Mach 0.85, AOA=0 deg, viscous case 
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Figure 18. Boundary Conditions definition. 
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a. The investigated domain (domain No. 1) was defined as a wall 

- Inviscid Adiabatic Wall (#1 in the table) for Euler solution. 

- Viscous Adiabatic Wall (#5 in the table). 



b. The two stings (Domains No. 2, 3) were defined as an Axis BCs 
(#15 - K around). 



No. 5) the body were defined as symmetry domains in the “z” 
direction, #13, ( i.e., the solution on +/- Z is identical), 
d. The outer shell of the block (Domain No. 6) was defined as 
Supersonic/ Subsonic inflow/outflow - #32. 

3. Reynolds (Re) Number Calculation 

In the case of viscous solutions, the Re Number must be calculated as follows: 
Assumptions: 

• Constant kinematics viscosity and equal to the value 



c. 



The two flat domains above (Domain No. 4) and below (Domain 



• Constant temp. T = 15 [C] = 

• Reference length - Body diameter L = 45.2 [cm] 



at 15° C 



v = 0.15 [cm 2 / sec] 

T = 15° [C] = 288 0 [Kelvin] 



Gas constant for air - 




Heat capacity for air - 



Nominal Mach Number - 



y = 1.4 
M„ = 0.85 



Re = — 



UL M„Lc M„LjyRT 



= 8.713*10 6 



V V V 



V 
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4. 



Results 



a. Mach Number Sensitivity: 

The nominal Mach number investigated was = 0.85 , which is 
approximately the critical Mach number. Since most of the research was done on 
transonic phenomena, this value was of interest. 

OVERFLOW is quite sensitive numerically at M „ = 1.0 . Previous 
researches by Biblarz and Priyono [Ref. 13] showed the difficulty of achieving good 
convergence at = 1.0. Hence, slightly different values were given; 0.98 or 1.02. 



The main reason for this sensitivity can be explained by the boundary 
condition implementation of the code. The code uses a method of characteristics to apply 
the boundary conditions. In general, as was described in Chapter n, five equations must 
be solved to get the 3 velocity components, the density and the internal energy. The 
eigenvalues of the equations are in the order of [Ref. 1 1]: 



A = 



U 

0 

0 

0 

0 



0 0 0 
U 0 0 

0 U 0 
0 0 U+c 

0 0 0 



0 

0 

0 

0 

U-c 



where U is the upstream velocity and c is the speed of sound under the specific 
conditions. U represents one of the boundary conditions defined to the block, the inflow 
condition. It can be seen that there is a singularity in the case where U = c and one of the 
eigenvalues becomes zero. 
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In this research, two sets of runs were done, for viscous and inviscid cases, 
both of which converged without any problems. Again, the reason for the ambiguity 
between this research results and Biblarz and Priyono’s research results might be due to 
grid sensitivities. 




Figure 19. CATM Cd pressure as a function of the Mach number. 

As can be seen from Fig. 19, the expected configuration of Cd as a 
function of the Mach number, where a maximum values are achieved slightly after M=1 
was not achieved. 
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This may be due to the following: 



• As was mentioned, the OVERFLOW code is sensitive to the transonic 
regime; hence, the results in this portion could be inaccurate. 

• Although the solution behaviors of convergence and stability were 
perfect, phenomena related to the grid generation may have caused the 
solution to be distorted. 

In the viscous case, no solution was achieved above Mach=1.5. 
OVERFLOW identified a negative density or pressure and stopped its iterations, and the 
“q.bomb” file was created for the last iteration. However, it can be seen that the pattern of 
the viscous case follows the pattern of the inviscid case. 

The solution problem probably originated in the grid configuration of the 
CATM. Since Mach=1.5 is far above the interested Mach regime, this phenomenon was 
further investigated. 

b. Solution Type - Euler vs. NS Solution 

Both the NS and Euler solutions converged smoothly through all the 
investigated Mach number regimes (0.5 - 1.8). The higher Mach numbers help the 
convergence to be more rapid. 

As mentioned, the difference between NS solutions and Euler was not 
dramatic; hence, most of this research was done with the Euler solutions, since they are 
much faster and consume less CPU. Figs. 20, 21 and 22 demonstrates the solutions at the 
nominal Mach number. 
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Cd pressure 



Cd pressure coefficient 




Figure 20. CATM alone- Cd pressure of Euler and NS solutions in M=0.85. 
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1 q- 6 Cd friction coefficient 




Figure 21. CATM alone, Cd friction of Euler and NS solutions in M=0.85. 
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Residuals - Delta Q 




Figure 22. CATM alone- residuals of Euler and NS solutions in M=0.85. 
c. Angle Of Attack Sensitivity 

As a technical remark, it should be mentioned that, since the orientation of 
the grid as developed using GRIDGEN was in a way that the symmetry plane was across 
the Z axis, the angle of attack and sideslip parameters had to be exchanged. The reason 
for it is that the angle of attack, a and sideslip angle, (3 are defined in OVERFLOW as: 



a = tan 



-i 



f \ 

w„ 



V«W 



( 20 ) , 



(3 = sin 1 




( 21 ) 
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where u„,v„,w m are the free stream components and U „ is the total free steam 
magnitude, equals to lu^.v^, w„| in the body frame [Ref. 10]. 

Hence, the definition of the AOA in the input file of OVERFLOW has to 
be under the sideslip angle (BETA) and vice versa, under angle of attack (ALFA), the 
value of the sliding angle has to be defined. 

Two main blocks were built; “jb_blockll.grd” and “jb_blockl3.grd”. 
The first was used mainly in the analysis, and the latter was a rotation of the first one by 
90 degrees to get a true orientation for a, (3 (see Appendix A for details). 

The sensitivity to angle of attack was investigated mainly on the CATM 
alone since the purpose was to analyze the influence of AOA separate from other 
parameters such as interference between the CATM and the missile. The sensitivity was 
checked in the NS viscous solution between 0 and 30 degrees. As expected, the drag 
increased rapidly with the AOA, see Fig. 23. 
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Figure 23. Cd of the CATM as a function of angle of attack. 

It should be mentioned that, in cases of convergence difficulties, 
increasing the angle of attack by several degrease (approximately 3-5 deg) helps to 
accelerate the solution convergence. This effect was published in Ref. 5 as well. The 
reason for it is because at zero angle of attack, the vertical velocity component is zero and 
the flow solver finds it more difficult to calculate the flow and the solution converges 
slowly. In higher angles of attack, other phenomena occurs such as flow separation and 
the convergence might slow again. This was not investigated in this research. 
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VI. MULTI-GRID CONFIGURATION 



A. GENERAL 

To investigate the combination of the CATM and the pylon, a multi-grid scheme 
was developed and run in various flow conditions and at relative locations between the 
two grids (see the Area-Rule in the following section). In addition, a wing grid was added 
to investigate the influence of wing existence above the CATM/Pylon system. 

It is common to use a peripheral grid (outer box) around the investigated grids in 
order to simulate the far field as closely as possible. This idea was applied in this research 
too; however, it was found that the CATM grid that was developed was broad enough so 
that the flow could fully conform to the free stream condition within the CATM grid 
itself, and an outer box was not needed. Comparison of the results with and without the 
outer box proved this; hence, the outer box was not used in further investigations. 

B. GRID GENERATION 

As a first attempt to analyze the flow around the CATM and the pylon, one grid 
was developed, consisting of the two configurations together. The intent was to prevent, 
at least for the first stages of the research, the necessity to deal with multiblock analysis 
which, as was described, is quite complicated. 

Mainly because of the differences between the length of the pylon and the CATM 
and the curvature lines of the pylon, no one satisfactory grid was found. The obstacles 
were “negatives cells”, recognized by OVERFLOW as preventing calculation of the flow 
by not updating these grid points in each iteration. An example of the combined block has 
been given in Fig. 3. In this block, as mentioned, the idea was to get rid of the negative 
cells by creating parallel grid lines with two identical domains: the CATM/Pylon 
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combination, placed apart from each other. Even this method failed, although logically it 
cannot create negative cells. 

The process of multi-grid generation was described previously. In this research, 
four separate grids were developed [Fig. 24]. 

• CATM grid - size 80 x 39 x 75 grid points. 

• Pylon grid - size 28 x 60 x 80 grid points. 

• Wing grid - size 80 x 30 x 39 grid points. 

• Outbox grid - size 100 x 80 x 60 grid points (not shown in Fig. 24). 

The outbox grid was not used exclusively in the process because, as was 
mentioned, the size of the CATM grid was broad enough that the flow got to its free- 
stream condition far from the boundaries of the block. 

One of the main issues of this research was the Area-Rule investigation. This 
subject will be elaborated upon in the next sections. In principal, the grid of the CATM 
was duplicated and shifted to get various cross-section area-distribution configurations. 

C. RESULTS 

The multi-grid cases were analyzed both in the CATM-Pylon-Wing configuration 
and the CATM-Pylon configuration only. As in the single grid, the research was done 
mainly in the transonic regime; hence, M=0.85 was chosen to be the main Mach number. 
The analysis was done both in Euler and NS solutions. Typical results are presented in 
Figs. 25 through 28: 
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Figure 24. Multi-grid configuration, CATM, pylon and wing. 
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Cd pressure 




Figure 25. Cd pressure development for the CATM, pylon and wing, inviscid flow. 
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Delta Q 




Figure 26. Residuals development for the CATM, pylon and wing, inviscid flow. 
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Figure 27. Cd pressure development for the CATM, pylon and wing, viscous flow. 
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Figure 28. Residuals development for the CATM, pylon and wing, viscous flow. 

It can be seen that, as was found in the single grid case, the difference between the 
viscous and inviscid solution is not a dramatic one. Typical Mach and pressure 
distributions are given in Figs. 29 through 32. 

In the following sections, the results are presented for the multi-grid analysis as a 
tool for the investigation of the Area-Rule. 
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Figure 29. Mach number distribution for a multi-grid case - NS solution. 
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Figure 30. Pressure distribution for a multi-grid case - NS solution. 
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Mach, invicid 





0.601388 



1.173005 



Figure 31. Mach number distribution for a multi-grid case - Euler solution. 
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Pressure, invicid 



/ 




Figure 32. Pressure distribution for a multi-grid case - Euler solution. 
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D. ORPHAN POINTS EXISTENCE SENSITIVITY 



1. General 

As mentioned previously, the appearance of ORPHAN POINTS is one of the 
main obstacles in the CFD process. As the number of grids increases, it becomes more 
difficult to get rid completely of the OPs. In some cases, a “no Orphan Points” situation 
has been achieved, but the combination of the grid was not physical, as will be elaborated 
as follows. 



2. Hole Definition and Orphan Points 

The definition of a hole and extrapolation regime was given in Chapter HI. In this 
research, when three grids merge together (CATM, Pylon and Wing), it was found that in 
order to get rid from all the Orphan Points, the limits of the grid of the CATM and the 
Wing to be calculated in the PEGSUS process (I,J,K INCLUDE parameter) can not 
include the whole range of grid points. 

3. The Physical Conflict 

The “physical” meaning of the hole definition is as follows: In the absence of 
Orphan Points, OVERFLOW calculates the flow over the CATM and the Wing as if the 
flow exists within a narrow portion in the Pylon (which makes the holes in the CATM 
and the Wing). This is a complete non-physical situation that can be identified by the fact 
that the flow pattern of the CATM and the Wing remain almost undisturbed on the 
symmetry plane of the blocks. 
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Attempts to get rid of all the Orphan Points and still preserve the physical 
condition of no flow through the Pylon failed. The attempts included, beside the 
techniques mentioned above, changes in the grids as well. 

4 . Sensitivity Study 

To investigate the physical necessity of the absence of Orphan Points, the 
PEGSUS input file was modified to prevent a “flow in the Pylon “ situation, while 
ignoring the existence of Orphan Points (which were found in a significant amount). It 
also was found that the aerodynamic coefficients almost do not change and the flow 
pattern becomes a “physical” one [Figs. 29 through 32], However, an optimization and 
sensitivity study should be done in any particular case to find the optimum situation of 
Orphan Points existence and to obtain a physical solution. 
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VII. THE AREA-RULE 



A. GENERAL 



Early in the 1950’s was found that the pressure drag of a specific configuration is 
highly dependent on the distribution of the cross-section area of the aerodynamic structure 
along the flow axis [Ref. 1]. The main pioneering work was done by Whitcomb and has 
been quoted since then in most of the literature. 



The benefits of the Area-Rule are significant mainly in the transonic regime. The 
concept can be demonstrated clearly by the missile and wings configuration (Fig. 33). To 
optimize the drag, the missile body was narrowed in the central portion and a “bottle 
neck” configuration was achieved, which compensated for the addition of the cross- 
section by the wings. 
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Figure 33. The use of the Area-Rule to reduce drag [Ref. 14]. 
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Another implementation of the rule was the increase of the cross-section area of 
the Boeing 747 passenger airplane at the front end. The central portion of the airplane, 
where the wings are mounted, was narrowed to achieve a smoother distribution of the 
cross-section area. 

However, the topic was mostly analyzed qualitatively or using parametric 
analysis. In this research an attempt was made to investigate and verify qualitatively the 
rule with CFD tools. 

A previous research at NPS, dealing with afterbody shape optimization [Ref. 15] 
found a strong dependence of the drag on the afterbody smoothness. The essence of the 
research is presented in Fig. 34 where it can be seen that, as the afterbody cross-section 
area develops in a smoother manner, the transonic drag decreases. 

0.045 1 1 —i — — ' 1 1 — — i 
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Free-Sfream Mach Number 



Figure 34. Drag Coefficient for Euler solution [Ref. 15], 
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Within all the written sources that were investigated in this research, the approach 
to the Area-Rule analysis was solely qualitative. It was obvious that the smoothness of the 
cross-section distributipn influences the drag directly and the smoother the distribution, 
the less the drag. Experimental results prove this theory. For example, in the technical 
note done by Lindsey in 1954 [Ref. 16] decreasing the height of the “bump” in the cross- 
section distribution plot created by the wing decreased the total drag of the body/wing 
combination. 

However, a quantitative or analytical analysis was not found. The main dilemma 
was the term “smoothness”. This subjective term can be interpreted differently by 
different individuals. Hence, in this research, with the use of powerful computers and the 
CFD process, the Area-Rule was modeled for quantitative analysis. 

B. THE AREA RULE IMPLEMENTATION 

1. General 

Originally, the main application of the Area-Rule was to combine slender bodies 
and wings to smooth the cross-section area distribution of the combination of the 
aerodynamic structures. In this research, the combination was different, the relative 
location of the CATM body and the missile was investigated to find whether changing in 
the cross-section distribution by translating the missile axial location backward and 
forward, relative to the pylon (Fig. 35), which can optimize the total drag of the 
configuration. 



$1 



Aligned Configuration 




Forward Configuration 




Backward Configuration 



Figure 35. Three basic models of CATM/pylon combination to investigate the Area- 
Rule. 

2. Method 

The influence of the cross-section area distribution was investigated by creating 
several combinations of pylon and CATM, of which each combination is different 
according to the location of the CATM relative to the pylon. 

A simple Matlab code was generated [Appendix 6], which loads the grid file as 
generated by GRIDGEN (after a few modifications such as deleting the “1” and the grid 
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size and completing the file to a full matrix by adding the appropriate zeros on the bottom 
of the file in order to be read by the code as a valid matrix). 

The code restriction is that the investigated shape must be grided geometry so that 
there will be constant “F domains. In the case of distribution along the X axis (or 
constant “J”, “K” in the case of distribution along the Y,Z, axis respectively). 

For this reason the pylon domain was hard to analyze, unlike the CATM itself 
where an exact distribution was found. 

The code calculates, for each cross-section, the cross-section area by integrating 
the area under the curve using the trapezoidal rule. 

It should be mentioned that each grid, even if it fulfills the restriction mentioned 
above, has to be carefully checked before running the code, since the structure of the grid 
file depends on the direction of the grid and the way the grid was created by GRJDGEN 
and modified by GRIDED. 

As a starting point of work, the following cases were analyzed: 

• Aligned case, i.e., the original relative pylon/CATM location. 

• Forward case, i.e., the CATM moved forward by 20” (-20”). 

• Backward case, i.e., the CATM moved backward by 20” (+20”). 

A distribution cross-section area is plotted in Fig. 36 and the three cases were run 
in OVERFLOW with and without a wing above. Although the smoothness of the cross- 
section distribution was improved in the backward case, there is a very sharp apex that 
causes a rapid change in the derivative. Hence a further optimization was done and the 
7.5 in.” was found to be the best in the aspect of “smoothness”. 
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Aligned case 




Forward case -20" 





Figure 36. Cross-section area distribution for CATM/pylon combination. 

C. INTERFERENCE DRAG 

Since we deal with a combination of a few aerodynamic configurations, we must 
consider the interference drag that occurs between the bodies in conjunction with the 
Area-Rule. 

One of the ways to control the interference drag is by “area ruling” the specific 
combination, i.e., by adjusting the cross-section area distribution to be as smooth as 
possible and, according to our research results, to control the slope of the afterbody. 
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Henfer and Bushnell give in Ref. 17 an example of a wing and nacelle, which is 
very similar to our case of the CATM and the pylon combination. 



The essence of the interference/area ruling can be presented in the charts shown in 

Fig. 37. 



Figure 37. Wing/nacelle example of drag reduction and interference drag [Ref. 17]. 

The results of this research described in Fig. 38, show a relatively high 
interference drag contribution. The drag of the CATM without the Pylon, i.e., without 
interference effects, is much lower then the drag of the various combinations of 
CATM/Pylon. 

D. AREA-RULE RESULTS 

1. No minimum value for drag was found for a specific cross-section 
distribution. 

2. The drag increases as the CATM location moves toward the backward of 
the pylon, and decreases as the CATM location moves forward, until it gets to the value 
of free-stream conditions. This can be explained by the lack of interference drag as the 
CATM moves far away from the pylon. 
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3. Since no special dependence was found for the cross-section area 
distribution, other “rules” were investigated, and it was found quite clearly that the drag 
of the CATM decreases as the back slope of the area distribution graph decreases. 

Using the code described above, several relative C ATM/Pylon location were 
investigated (aligned case, +20”, -20”, -10”. -7.5” and -5”). For each case, the area cross- 
section distribution was evaluated and, using linear regression, the slope of the back part 
of the graph was evaluated. The results are presented in Fig. 38 and include the reference 
of the free-flight case. 



x 10 3 




Figure 38. Cd as a function of the average slope of the afterbody cross-section area 
distribution. 
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VIII. THE EQUIVALENCE RULE 



A. GENERAL 

The Area-Rule, can be considered an extension of another significant rule, the 
Equivalence-Rule [Ref. 14]. 

In essence, the equivalence rule is a tool for treating relatively complicated 
aerodynamic configurations by comparing the investigated configuration to a simpler one 
that has a more trivial solution, for example, a slender body of revolution. 

Particularly in the transonic regime where the emphasis of this research lies, one 
can reduce the linearized, general compressible flow equation [Ref. 14]: 

(l-M ; )®„ +<!>„+ ®„=0 (22) 



to the simple Laplace equation: 



4> yy +<D Z2 =0 (23) 

and apply the perturbation theory on the simple model. 

In subsonic and transonic flows, the far-field flow structure depends on the body 
shape and volume, but in the transonic regime the far field depends on the volume only. 
On the other hand, in transonic flows close to Mach 1, the waves move out in a 
perpendicular direction to the flow. Hence we can deduce that, for the transonic regime, 
the far-field flow depends on some configuration, not necessarily the actual, investigated 
one, that has the same cross-section area distribution as the original geometry. 
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Naturally, we would prefer to choose a simple body of revolution that has the 
same cross-section area distribution and is easy to deal with in a computational point of 
view. 



B. SMALL PERTURBATION THEORY 

The assumption that the equivalence theory is based on the small perturbation 
theory [Ref. 14] leads to the relatively simple irrotational flow. Its velocity would satisfy 
the gradient of the potential 4>(x,y,z) which can be represented, in a zero angle of attack 
as: 



<E>(x, y , z) = U„ (x) + 0(x, y, z) (24) 



where <t>(x,y, z) satisfies the following differential equation: 

(1-M„ 2 )4»„+4» w + 0« =Mj[(y + l)/^]0,4»„ (25) 

C. THE AREA-RULE WITH THE EQUIVALENCE BACKGROUND 

Heinrich [Ref. 14] discussed the equivalence mle in detail and derived the exact 
expressions for the solutions of the basic perturbations equation for the transonic regime, 
where: 



= (26) 
which is similar to the heat equation and is its solution. 
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The solution for the transonic regime will be: 



i = ™ h !£ + I|'M 4 (27) 

U 4n 4nx 4n ™ x — £, 



and for the equivalent body: 



u S"(x), 5 

— = In — 

U 4n n 



(28) 



The smoothness of the cross-section area can be well represented by the behavior 
of the second derivative of the cross-section area function, S(x ) = n y b . A sudden 
change in the cross-section area function will increase abruptly the value of S" (x) and 
hence increase both the drag of the equivalent and the original body. 



The subtraction in the last equation represents the contribution of the Cp which 
originated in the cross-section area distribution along the body axis. 



D. SLENDER BODY THEORY 



The Area-Rule can be most useful when it is applied on a simple configuration as 
a slender body. Finding a slender body with a cross-section distribution that is similar to 
the investigated configuration, can make the transonic drag analysis much simpler. 

The slender body theory [Ref. 14] states that the flow potential around a slender 
body satisfies the following: 



O = <h 2 + g(x) (29) 
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where <b 2 is the solution for the 2D degenerated equation (Laplace’s form) and g(x) is an 
addition function that depends on the free stream condition and the cross-section area 
distribution . 
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Figure 39. Illustration of the transonic equivalence rule [Ref. 14]. 



The pressure distribution on a slender body of revolution is as follows [Ref. 14 
and 17]: 



C P = P P ~ =-—<£, — -(<b 2 +<b 2 ) = 

F r\ IT 2 / O TT U ^ y ' 



P-UJ/2 U 



tm: 



i+ 



J -Ml 



t (U,+ b) 2 +v 2 + 
U„ 



w 



2 \ 



Y-l 



(30) 



-1 



After further simplification, assuming the square of the velocities in the y, z are 
much smaller then in the x direction, and after applying binomial expansion, we get: 
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Obviously, in the transonic regime, the x term can be neglected and the following 
expression for the pressure coefficient can be obtained: 



C P 

1 Transonic 




2 2 

v + 



U 



(32) 



From here, the total forces on a slender body can be found and, using the 
combination of the area rule and the equivalence rule, the forces on the investigated 
configuration can be evaluated. 

E. THE TRANSONIC SIMILARITY RULE 

1 . The powerful ability of the similarity rule can be found specifically in the 
transonic regime where the flow over the body reaches the M = 1 limit. The main obstacle 
is to measure and calculate the aerodynamic coefficients. The similarity rule makes it 
possible to use data recorded for subsonic and transonic regimes, and interpolate it for the 
problematic transonic regime. 

2. In the transonic regime, when the first shocks are being built, one can 
assume that the gap from M = 1 is identical before and after the shock [Ref. 14], If we 
denote I, II as the conditions before and after the shock, we will get: 

M,— 1 = 1 — M n (33) 
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Figure 40. Transonic wave drag [Ref. 17]. 



The pressure coefficient derivative with respect to the free-stream Mach number, 
can be found from the expression for Cp. Assuming the slender body and small 
perturbations theory, we get: 
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(for slender body) 



(34) 



The derivation for the drag and lift coefficient is clear, and will not be presented 
here, although the final implementations can be seen in Fig. 41 [Ref. 14]. For a wedge of 
thickness h, the drag coefficient at M=1 can be found for various wedge angles and 
generally, as was found by Spreiter [Ref. 18], as a function of the similarity variable %. 
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(a> 




Figure 41 . Cd as a function of free stream Mach number - and similarity variable-% 
[Ref. 14]. 

F. COMPARISON TO OUR RESULTS 

In addition to the smoothness influence as was found by Whitcomb [Ref. 1 ] and 
reported by others such as Heinrich [Ref. 14], it was found in this research that the cross- 
section area distribution influences the drag behavior, not only according to the 
smoothness of it’s development, but according to the slope of its development along the 
afterbody. 
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IX. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 



A. GENERAL SUMMARY 

1. A CFD solver was applied to a multi-grid configuration and the precise 
procedure has been documented in detail so as to assist further projects. 

2. CFD process was found to be very sensitive to the grid configuration. The 
grid architecture, including grid line distribution, blocks structure, and grids merging, 
significantly affects each stage of the process, particularly OVERFLOW. 

3. A primary dependence of the drag on afterbody average cross-section area 
distribution slope was found which can guide aerodynamic shape development for 
minimum drag. This addition to the traditional Area-Rule can be a powerful tool to 
accelerate development time. 

B. CONCLUSIONS 

1 . The minimum drag of the CATM mounted on a pylon was found when the 
CATM is shifted 20” toward flight direction (-20” configuration) and its value is 
0.003766. It can be seen (Fig. 38) that the drag coefficient reduces as the CATM gets 
further from the pylon toward the free flight condition. 

2. The drag behavior of the CATM/pylon combination can be explained 
mainly due to the interference drag caused by the CATM/pylon relative location. 

3. Area Rule was found to be dependent on the afterbody cross-section area 
distribution average slope, and not on the overall smoothness of the cross-section 
distribution of the CATM/Pylon configuration. The average afterbody cross-section area 
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distribution correlates well with the drag of the different relative location of the 
CATM/pylon 

4. The results of Euler solutions rather than NS solutions are, in most cases, 
quite adequate and save the computational efforts. 

5. Using CFD, Similarity-Rule can was verified only qualitatively, further 
work has to be done to analyze the Similarity-Rule quantitatively as was done with the 
Area-Rule. 

6. Orphan Points can be found in the multi-grid approach, but are tolerable as 
long as their existence does not affect the physical meaning of the solution. In any case, 
the preferred solution contains no Orphan Points, but, achieving this is not always worth 
the effort. 

7. The use of the NPS Cray services as a main operating system causes 
significantly delays and increases time response for the results. The optimum and CPU 
time should be developed economically in order to minimize negative effects, either by 
using the correct “*.nqs” file or by carefully selecting the critical parameters to achieve 
rapid convergence. 

C. RECOMMENDATIONS 

1. To simplify the work in modeling the CATM configuration, and make the 
results more accurate, a current geometrical data base should be obtained from the 
contractor (TI) that can be read by the grid generation code in NPS, such as GRIDGEN. 

2. The best location to minimize drag is where the CATM is as far forward as 
possible on the pylon. This result affects structural considerations on the one hand, and 
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antenna locations (GPS on the-front end of the CATM and communication on the rear- 
end) on the other. A CATM placed too far forward will cause difficulties for the guided 
aircraft because of antenna blocking. 

3. The basic study of any multi-grid case configuration should start by 
studying each of the grids separately to assure valid grids that can produce physical 

solutions. 

4. When dealing with multi-grid cases, an outer box should be used unless 
the original, investigated blocks are as broad that the effects from the body do not 
propagate out to the block boundaries. Therefore, verifying the characteristics of the 
separated grids, should be done in case an extra block has to be added that will result in 
significant increases in the computational needs and complexity. 

5. Since OVERFLOW also runs on SGI machines, it is recommended to run 
the initial, non-CPU consuming cases on the SGI machines and to check them directly by 
FAST (i.e., “q.save” generated by OVERFLOW, which runs on the SGI can be read 
directly by FAST as a “solution”, “unformatted” file, on the other hand, in the Cray case, 
the “q.save” file has to be converted to a SGI readable format using a “wrq” conversion 
file). 



6. The complexity of the CFD process can be eased by making the codes 
more user friendly. Recommendations for such changes are given in Appendix E. 

7. Since grid generation is one of the most difficult phases in the CFD 
process, a reliable grid generation tool is required. GRIDGEN is an excellent tool, and it 
is recommended that the most current version be installed including all the technical 
support that accompanies it. 
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8. As mentioned in chapter HL-F, the GRIDED version currently installed is 
not the latest version which has more useful applications. It is recommended to get the 
updated version, which can help to create simple blocks easily by modifying 2D grids. 

9. A further study of a full JSOW including wings and fins is recommended 
in order to analyze the aerodynamic characteristics of the operational weapon system to 
be used in other applications such as guidance simulation. 

10. Some work should be done to analyze the store-separation behavior of the 
CATM. Using two codes that were written by Nielson Engineering, SUBSAL and 
SUBSTR, one can investigate the store trajectory after modeling the store as an 
equivalent axisymmetric object. The codes can be provided by NAWC China Lake and 
are ready to use on the SGI machines. This process might complete (partially) the picture 
by analyzing the CATM separation from the F- 1 8 as a model carrier. 



98 



LIST OF REFERENCES 



1. Whitcomb, T.R., NASA Report 1273 (1952); also, Aviation Week and Space 
Technology, Sept. 19, 1955 and IAS paper 601 (1956). 

2. “US authorized two further JSOW versions”, IDR International Defense Review - 
Jane’s Magazine, July 1995. 

3. Fulghum, David A., “New Wartime Roles Foreseen For JSOW”, Aviation Week 
and Space Technology, February 28, 1994. 

4. John, D., “NAVY Primed for JSOW critical Design Review”, Aviation Week and 
Space Technology, February 27, 1995. 

5. Kern, Steven B. and Bruner, Christopher W., “External Carriage Analysis of a 
Generic Finned - Store on the F-16 using USM3D”, AIAA-96-2456-CP. 

6. Rom, Josef , “ High Angle of Attack Aerodynamics” , Springer - Verlag,1991. 

7. Barton, Bret S., “Application of Multi- Block CFD Techniques to a Missile 
Geometry”, Master of Science in Aeronautical Engineering, Naval Postgraduate 
School, December 1995. 

8. Steinbrenner, John P. and Chawner, John R., “GRIDGEN Release Notes: Version 
9.6 (First Production Release)”, MDA Engineering, Inc., October 18, 1996. 

9. Suhs, N. E. and Tramel, R. W., “ Pegsus 4.0 User’s Manual’, Calspan 
Corporation/ AEDC Operations, June 1991. 



99 



10. Buning, Peter G., NACA Ames, “ Overflow User’s ManuaF, Version 1.6aw 12, 
July 1995. 

1 1. Merkle, Charles L., “Computational Fluid Dynamics of In viscid and High 
Reynolds Number Flow”, AA4507 notes. Fall September, 1990. 

12. Scarry, Michael T., “JSOW CATM Conceptual Weight and Aircraft Design ”, 
Master of Science in Aeronautical Engineering, Naval Postgraduate School, 
September, 1996. 

13. Biblarz, O. and Priyono E., “Transonic Pressure Drag Coefficient for 
Asymmetric Bodies”, Proceedings, ICAS-94-2.5.2 19th Congress ofICAS, 
Anaheim, CA, 1994. 

14. Heinrich, J. Ramm, “ Fluid Dynamics for the Study of Transonic Flow ”, New 
York - Oxford, Oxford University, 1990. 

15. Priyono, Eddy, “An Investigation of the transonic pressure drag coefficient for 
axi-symmetric bodies”. Master of Science in Aeronautical Engineering, Naval 
Postgraduate School, March 1994. 

16. Lindsey, W.F., “The Flow Past An Unswept And A Swept Wing Body 
Combination And Their Equivalent Body Of Revolution At Mach Numbers Near 
1.0”, National Advisory Committee For Aeronautics, Technical Note 3703, 1954. 

17. Hefner, Jerry N. and Bushnell, Dennis M., “An Overview of Concepts for Aircraft 
Drag Reduction”, NASA Langgley Research Center, Hampton, Virginia, 
AGARD-R 654, 1977. 



100 



18. Spreiter, John R. and Stahara, Stephan S., “Aerodynamics of Slender Bodies and 
Wing-Body Combinations at = 1 ”, Nielsen Engineering & Research Inc., 
Mountain View , Calif., AIAA V 1 9, No. 9 Sep 1971. 

19. Fan, Yue Sang, “An Investigation of The Transonic Drag Coefficient for Axi - 
Symmetric Bodies”, Master of Science in Aeronautical Engineering, Naval 
Postgraduate School, March 1995. 



101 









102 



APPENDIX A. GRID FILES DESCRIPTION. 



A. GENERAL 

The evolution of this research is reflected by the grids created throughout the 
work. Since the major obstacles were found in the grid generation and caused by grid 
structures, we found it very important to keep a detailed description of the grids and their 
related phenomena, although most of them probably will be erased after the completion of 
this research. 

Generally, each grid was modified after it was created to achieve a uniform 
configuration. The modification was done by GRIDED and the modified grid was 
distinguished from the original one. The original grid has a “.grd” suffix, the same grid 
without the “1” on the top of the file (which indicates for GRIDGEN that it is a one-block 
grid and has to be eliminated in order to be read by GRIDED), has the “_m.grd” suffix. 
The same grid after being modified by GRIDED has, in addition a suffix indicating the 
changes done by GRIDED, e.g., “_4_100.grd” means the grid changed J, L indices 
(option “4” in GRIDED) and changed the direction of J index. Finally “uf ’ suffix indicates 
that the file is unformatted (and can be used directly for OVERFLOW needs). 

The following description includes only the basic configurations. 
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B. THE GRIDS 



I. Test case grids 

a. testgrd - 80x39x75 

This grid was generated to study the basic concept of OVERFLOW after 
attempts to get a valid solution for the CATM itself failed because of consistent 
unexplainable negative values of Cd for subsonic flow and positive for supersonic, as the 
CATM grid, it is a 80x39x75 grid, symmetric on Z axis. The purpose was to create the 
simplest grid that won’t make any special difficulties due to asymmetry, non slender etc. 
No good solutions were achieved using this grid, probably because of two reasons. 

• Although the grid was symmetric, it was not a body of revolution; 
therefore relatively sharp comers were formed, including on the plane of 
symmetry. The thought was that comers were not handled by denser 
grid lines might cause divergence of the solution. 

• The other suspected reason is the distribution of the grid lines from the 
surface outward, although a geometric distribution function was applied 
(AS Initial -0, final 20). The grid was not dense enough, or the AS Initial 
-0 input created difficulties. 

b. testl.grd - 80x 7 7x 75 

One reason for the failure of the test.grd file may be that the boundary 
conditions definition was incorrect, especially the symmetric BCs. In order to change the 
BCs radically and check the theorem, a full size grid was created (not a half 
configuration). 
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The grid was not used, since, it creates negative volumes in OVERFLOW 
with no ability to overcome them. 

c. testjdensedgrd - 80x39x75 

The same as test.grd but with a different grid distribution , AS Initial -0, 
final 50. It appeared that even this distribution did not help and the negative Cd continued 
to appear. 



d BasebalLgrd - 80x39x75 

Since the test.grd grid was not a body of revolution, a similar grid was 
created, but this one was a body of revolution. The change was done in order to 
investigate the “kinks” in the grid, on the symmetry planes in the solution. No 
improvement was found, and the same negative Cd occurred, which indicates that these 
“kinks” have no significant influence on the solution. 

e. bball_small.grd - 80x39x75 

Looking at a few sample cases, such as the Buning sample case and the 
Barton thesis case, one can see that the width of the blocks is narrow (relative to our 
blocks, which were made very broad in order to prevent the influence of the block 
boundaries and reflections). The narrow blocks were used because the cases were 
multiblock, and on the top of these blocks was an “outerbox” to calculate the far field 
conditions. However, investigating the grids as a monoblock case, one can find a valid 
solution. Therefore, a very narrow block was created while the grid size remained the 
same. 
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No good solution was found, but the phenomena of negative cells 
disappeared, and a stable positive value was found, with the restriction of high residuals 
values (approximately 10E-01 to 10E-03). 

f. bb_dense.grd - 80x39x75 

After the origin of the negative Cd was found, a new grid was created, 
identical to the “baseball.grd” but with a different grid distribution. A very dense grid was 
created near the body surface. To be more specific, a geometrical distribution function was 
chosen with a AS initial and final of 0.05 and 100, respectively (a distribution of 0.01 - 150 
failed due to negative cells in OVERFLOW). 

Negative Cd was not found, but the lack of convergence was found. 

g. Baseball_90.grd - 80x39x75 

To achieve a maximum computability between the test case and the Buning 
test case, the baseball.grd grid was rotated by 90 degrees, so that the symmetry plane was 
“Y” instead of “Z”. The solution, as expected, was the same. 

tu bb_dense_ang.gr d - 80x39x75 

In a few projects [Ref. 7 and Ref. 19] and in Buning sample case as well, 
the bodies investigated were not a half body but slightly more, i.e., approximately 200 
degrees were investigated rather then 180 degrees as expected. Therefor a new grid was 
created out of the basic “baseball_90.grd” grid but in 200 degrees. No different 
phenomena were found. 
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The reason for the addition of the angle was not clear, especially since the 
BCs on the inclined surfaces is a symmetry in the “Z” direction and, once the surface is 
inclined, the symmetry is not purely in the “Z” axis. 

2. JSOW/CATM Grids. 

a. Jbjblockd.grd - 80 x 39 x 85 

This is the basic JSOW block that was created. It is a very narrow block. 
The diameter of the outer shell is approximately twice the length of the missile itself. 
Hence, phenomenon of shock reflection occurred and can be seen in Fig. Al. The reflected 
shocks hit the body near its rear end and create a “ring” that can be seen in Fig. Al 

b. Jb_block5.grd - 80 x 39 x 85 

Due to the shock reflection phenomenon, the block was significantly 
increased and, indeed, the reflections disappeared. 
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Figure Al. Reflection of shock waves due to too small block size. 
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c. Jb_block8.grd - 80 x 39 x 75 

This is the first short version of the CATM after it was modified and its 
weight was decreased to below 300 lb [Ref. 12], The grid size decreased a little in the “K” 
direction. No special changes were made to the grid distribution. 

d. Jb_block9.grd - 80 x 39 x 75 

This is the same grid as the jb_block8.grd, but with a wider block, created 
to get rid of the “ring” that appears at the rear portion of the missile due to a reflection 
from the block boundary [Fig. Al]. The “ring” phenomena disappeared, but, the solution 
gave negative Cd for the subsonic regime (Mach number of 0.5 to 0.8) and a positive Cd 
in the supersonic regime (a Mach number of 1.2). 

e. Jbjblockl O.grd - 80 x 39 x 75 

This is a very wide block (radius to body length ratio of about 1500), was 
created since it was suspected that the previous block was not broad enough 
(approximately 3 to 4 times the body length) and there still is an influence of the outer 
shell boundaiy on the solution. It was found that it is not true, and the block does not have 
to be so wide, especially in a multiblock grid. The block gave consistently negative Cd. 

f. Jb blockll.grd - 80 x 39 x 75 

This is a block that was created after it was found that the dense grid near 
the surface is more important than was initially thought. The distribution function was set 
to “geometric” with an initial and final spacing of 0.05, 150. The negative Cd phenomena 
did not occur. 
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The grid was modified again and an extra plane was added to the J domains 
on the other side of the symmetric plane (plane Z=0). The addition of the plane was done 
by using the GRIDED feature to add another domain. The purpose of the extra plane is to 
enable OVERFLOW to imply the BCs on the symmetry plane in so that there will be 
“sources “ of information from both sides of the plane to get by averaging the information 
from the middle (symmetry) plane. 

k JbblockU.grd -80x39x 75 

Before the “extra” plane was added to the previous block 
(Jbblockl l.grd), an attempt was made to run a case with a block that was extended to 
more then 180 degrees. Obviously it was a failure, since only one surface had to be added 
(as described above) and, in this block, more then one grid was added. 

i. Jb_blockl3.grd - 80 x 39 x 75 

Since the Jb blockl l.grd block was oriented with the Z axis as the 
symmetry axis, the definition of a, (5 was disturbed; hence, an identical block was created 
that was rotated by 90 degrees relative to Jb blockl l .grd in order to verify the definition 
of the directions of a,p. 

Since all the research was done using the Jb_blockl l.grd block, and the 
orientation of the pylon and wing blocks was similar to the Jb_blockl l.grd, the 
Jb_blockl3.grd block was used only to verify the directions of a,p. 

3. Pylon Grids 

The Pylon grid originally was generated as the JSOW/CATM from a set of 
drawings. Since the main reason for using the Pylon is to investigate the Area Rule rather 
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than to analyze its aerodynamic behavior, it is not critical that the geometry be identical to 
the reality. 

The shape of the Pylon was given as contour lines only; the details of the geometry 
was unknown and assumptions had to be made. Many Pylon grids were created; however, 
unlike the JSOW/CATM case, this block was not a problematic one and through most of 
the research only one block ( “p_densed.grd”) was used. The block was modified in 
several cases by applying the elliptical solver in GRIDGEN to prevent Negative Volumes 
(p_dense_elip5.grd - 5 iterations of elliptical solver and p_dense_elip8.grd - 8 iterations of 
elliptical solver). 

In the research phase, where Orphan Points were found when a combination 
between the CATM and the Pylon grids was done, a new block was created - 
(“P_den_ang.grd”) which had a different block structure than the original. If the original 
block (Fig. 24) was built as half of cylinder, then in the “angle” block, the upper and lower 
domains were tilted 45 degrees to “open” the block. The reason for doing so was to create 
grid lines that will spread at an angle and will be more parallel to the grid lines of the 
CATM. Thus the grid cells of the Pylon far from the Pylon surface will be similar in size to 
the CATM grid cells in the same area and, the information from the Pylon to the CATM 
would be propagate in a better way. The modification helped slightly to get rid of part of 
the OPs; however, the block was not used further because other methods, as described 
previously, were applied. 

4. Wing Grid 

The wing grid was created to investigate the sensitivity of the flow pattern around 
the CATM and under a wing. Since no accurate information about the geometry of a F- 18 
wing was available, a general wing grid was created, “w_blockl_m.grd. The block was 



111 



similar to the CATM block, an “0” shape block with two stings and two symmetry 
domains. 

5. Outer Box 

As was mentioned previously, the outer box was created but not used much since 
the CATM block was broad enough that the influence of the body did not propagate to the 
block boundaries. 

The block grid, “outbox_bl.grd”,was a simple cubic shape block, in order to save 
grid points and CPU, the distribution of the grid points was set so that the grid was dense 
around the CATM / Pylon and the grid became coarser toward the outer boundary of the 
outer box. 
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APPENDIX B. FILES OF INTEREST IN THE CFD PROCESS 



A. GENERAL 

This Appendix includes several files: executables and input files that are being 
used in the CFD process and can be used as an example for future research. 

B. EXECUTABLE FILES: 

• xa2b.f 

• xb2a.f 

• rmg2peg.f 

• merge41.f 

C. INPUT FILES 

• OVERFLOW input file sample. 

• PEGSUS input file sample. 

D. INFORMATION FILES 

• ja_sum.out 
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APPENDIX C. CFD PROCEDURE FLOW CHART 



(Based on the flow chart in Ref. 7) 
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APPENDIX D. SCRIPT FILES 



A. GENERAL 

This appendix includes a sample of script files created during this research to 
make the CFD process more efficient and “user friendly”. The script files are important 
because of the various files and processes that must be run many times in the same 
manner. Using a script file assures that the process will be done in the right order. 

B. THE FILES 

1. scriptl. A script file used to run RMG2PEG, PEGSUS, MERGE41 and 
OVERFLOW. Originally the file was created by Buning and can be found in the sample 
cases of OVERFLOW. It was also used by Barton [Ref. 7]. 

In this research the file was modified to handle more codes and, hence make the 
process more efficient. 

2. results. A script file created to run the Matlab codes READFOMO.M, 
READRESED.M and to plot the Cd and residuals. 
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APPENDIX E. OVERFLOW USER’S INTERFACE IMPROVEMENT 



SUGGESTIONS 



A. GENERAL 

OVERFLOW is a code that is still under development. Few recommendations are 
given to the user in advance to improve the “user friendly” level of the code especially 
concerning error warnings. If these modifications can be made to the code, future users 
can save a significant amount of time and frustration in their research. 

B. RECOMENDATION FOR IMPROVEMENTS. 

1. Symmetric Plane Boundary Conditions. 

BCs number 11, 12 and 13 imply symmetry BCs for the X,Y and Z directions 
respectively [Table 3.3, Ref. 10]. To use these conditions one must add an extra plane so 
that the symmetry plane will be in between two grid domains and its information will be 
calculated from these two domains. This action, as detailed earlier is done by GRIDED. 
However, if this extra plane is not added, a valid (“pseudo”) solution is given, which is 
obviously not physical. 

An error message should be given in the OVERFLOW output file if the code 
recognizes a definition of the symmetry plane without adding an extra plane. 

2. Distinguishing Between Equal Grid Size Grids. 

Two grids of the same size can provide the same solution “q.save” file, no matter 
what their shape. The code has to recognize a valid grid and solution pair and prevent 
mixing a grid and a solution from another grid. 
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This improvement is important when dealing with a significant number of grids 
that are only slightly different one from the other (e.g., in the number of Elliptical Solver 
iterations implied on the grid). 

3. Mixing Between NS And Euler Solutions. 

The difference between a NS and Euler solution is physically dramatic. However, 
to flip between the two in the OVERFLOW input file is a very simple command. One can 
easily change the input file to an NS or Euler solution by changing the VISCJ, VISCK 
and VISCL parameters in the VISINP group. 

To get a valid solution, the BCs must be changed as well (e.g., for a wall, from an 
inviscid adiabatic wall to a viscous adiabatic wall). If not - a solution will be resolved, but 
it will not be a physical one; hence, an error message should be raised if there is a 
mismatch between the VISINP data and the BCs type. 
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APPENDIX F. MATLAB CODE GENERATES AREA CROSS SECTION 



DISTRIBUTION 
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APPENDIX G. FEW RECOMMENDATIONS FOR THE CRAY USERS. 



Most of the jobs, especially the heavy multiblock cases, were run on the Cray 
computer in the Computer Center of NPS. The interface format of the user with the Cray 
is by job submitting. 

The current Cray model is Y-MPEL89, and consists of 8 CPU’s sharing 256 
Mwords of memory. 

As for this research, the Cray (Jedi) was able to deal with only ten jobs 
simultaneously. In any case, a user cannot run more then one job at a time. The jobs are 
sent by “nqs” files [App. H], which define the job type, commands, and priority. The 
only approved priority that can be used by the students are the “reg” or “econ” types. The 
“prem” type is forbidden to be used. Because of the this, jobs can wait for a long time to 
be executed after being submitted by the user. 

To avoid too much waiting time when the user identifies a “busy” Cray (by typing 
the “nqstat” command in the Cray), the runs can be divided into relatively low numbers of 
iterations, and the “restart” feature of OVERFLOW can be used to define low amount of 
memory required for the each job. 

The job’s priority decision is done, not according to “first submitted, first 
executed” but, according to the amount of time and memory required as defined in the 
nqs file. Long, “heavy” jobs wait in the queue a longer time. 

A practical method to have an idea about the size of a job, how much time it runs, 
how much CPU it consumes and more is through “ja_sum.out” file. This file contains 
some useful information that might help one to plan his jobs. An example is given in 
App. B and the command lines to get the file are given within the “nqs” file in App. H. 
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APPENDIX H. NQS FILE SAMPLE 



A. GENERAL 

Since the interaction with the Cray is through batch files and not an “on line” 
process as is usually performed, an “nqs” file has to be created, in which the procedure of 
the commands is given. 

The benefits of the nqs file were not fully used in this research. This file can be 
used to create a complicated procedure, including running jobs, transferring files, and 
changing directories. 
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the block (e.g., M(+20”) ). 
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Reynolds number and Temperature at infinity values 
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APPENDIX J. A PRINCIPAL PROCESS OF CFD ANALYSIS 



A. General 

A concise flow chart describing the process of CFD is given in Ref. 7. The chart 
includes most of the phases in the process; however, it was found that a detailed 
description of the different stages including some emphasis on particular important points 
can be very helpful for the user. 

Therefore, it is recommended to use the flow chart, given in Appendix C parallel 
to the literal description. 

The description will concentrate on the multiblock case as was done in this 
research (i.e., the tools that we used, e.g., GRIDGEN and FAST). From this general case, 
it is obvious how to handle the single block case. 

B. The Process 

1. Generate the grids (using GRIDGEN or any other grid generation tool). 

2. Edit the grids and remove the “1” from the top left comer of the file. This 
number describes the number of blocks included in each file, and has to be removed in 
order to ran GRIDED and “xa2b” on the file. It is recommended to rename the “1 less” 
file in order to prevent long editing times in case the need to read the grid in GRIDGEN 
arises again. 

3. Visualize the grids using FAST; write down the directions and orientation 
of I, J, K in order to plan the changes that have to be done on the grid. 
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4. Decide on a particular grid direction and orientation that will be consistent 
through the research. That way, only minor changes will have to be made to the 
OVERFLOW, PEGSUS input files, and quite a few sources for mistakes will be 
eliminated. In our case, for example, through all the research, the missile surface was 
K=l, the front sting was 1=1 and the lower symmetry plane was J=l. 

5. Run GRIDED on the file and change the orientations and directions of the 
grid as needed. In addition, add an extra plane in the symmetry planes. 

6. Convert the file created by GRIDED to a binary format using “xa2b” code. 
Pay attention that the file has to be compiled differently if the process will continue to be 
done on the Cray or locally on the SGI machines. 

7. In this stage, it is recommended to run each grid separately on a simple and 
short (approximately 50 to 100 iterations) OVERFLOW case to check, in principle the 
convergence of the solution and absence of negative cells. If a negative cell error were 
found in the OVERFLOW output file, before going back to the long process of grid 
generation and modification, it is recommended to run GRIDED again on the original file 
(i.e., do steps “b” and so on again), and/or compile the conversion file “xa2b” and run it 
again. If this does not help, a modified grid has to be generated with different grid line 
spacing (distribution) and/or with an elliptical solver implied on the domains. It was 
found that negative cells appear sometimes without any ’’physical” reason. 

8. Prepare the grids to run PEGSUS by applying the RMG2PEG code on the 
grids. It is very convenient and saves time if a script file is generated to run the process of 
RMG2PEG -> PEGSUS -> MERGE41 -> OVERFLOW (Appendix D). Run the script 
file for a small amount of OVERFLOW iterations to verify the whole process. Check 
Orphan Points in the PEGSUS output file. It must be mentioned that OPs will not prevent 
a valid OVERFLOW solution. As mentioned before, OPs is a point that is not being 
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updated; however, it is recommended to get rid of all the OPs to prevent “holes” in the 
solution. In the case of OPs, change the PEGSUS input file according to the 
mechanization’s detailed earlier. 

9. Once a good PEGSUS input file has been achieved, OVERFLOW can be 
run for a significant number of iterations. The physical solution can be analyzed starting 
from residual values of approximately 10E-5 (AQ); however, the other way around is not 
true, these residual values do not imply a good solution. A physical value of aerodynamic 
coefficients and a reasonable flow pattern has to be found. 

10. The solution can now be visualized using FAST, in the case of. a Cray run, 
the “Q” file, “q.save”, has to be converted to a SGI format using the “wrq” file, and the 
“grid.in” file has to be converted to a SGI format using the “wrgr” file. Notice that the file 
created by “wrgr”, “gr.sgi”, includes the blanking points as well as the whole other grid 
points. Therefor, if one wants to view the blanking scheme, that file has to be used in 
FAST. In the case of a local, SGI run, the solution “q.save” can be read directly as a 
solution/unformatted file. 

1 1 . More OVERFLOW iterations or modifications of the grids can be done on 
the base of the successful run. 
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